2017/9/07
DBDA2E-utilities.R # 便利関数 Stan-Ydich-XnomSsubj-MbinomBetaOmegaKappa-Power.R # 第13章タッチセラピーの例数設計 Stan-Ydich-XnomSsubj-MbinomBetaOmegaKappa.R # タッチセラピーモデルのベイズ推定コード
著者の自作関数を読み解くとこが結構疲れる
想定はStan中級者 (階層モデルはいけるぜ!)
(Just how many times must I show her I care,)
(Until she believes that I'll always be there?)
(Well, while she denies that my value's enough,)
(I'll have to rely on the power of love)
と言いつつ…
愛の力=\(検定力\) 僕の価値=\(目標とするパラメータ値\) あの子=\(査読者\)
| 帰無仮説が真(差なし) | 帰無仮説が偽 (差あり) | |
|---|---|---|
| 検定結果偽(差なし) | 1-α | β |
| 検定結果真(差あり) | α | 1-β |
Note: α= Type 1 error, β= Type 2 error
例えば、歪んたコインで表が出る確率がチャンスレベル(50%)の周辺から外れているか
事後分布の幅の狭さ
予測値の尤もらしさ、予測区間
空値周辺の実践的に等価な範囲(ROPE)が事後分布の95%HDIを含まないことを示す
事後分布の95%HDIが特定の最大幅より狭いことを示す
事後分布の95%HDIを含む予測値のROPEを示す
事後分布の95%HDIの下限が、ROPEの上限以下
事後分布の95%HDIの上限が、ROPEの下限以下
事後分布の95%HDIが、ROPE[-0.7,0.7]以内
事後分布の95%HDIが特定の最大幅より狭いことを示す
2000回投げて,65%が表になる
#仮説パラメタ分布のモード (ω)と集中度 (κ)を指定 kappa= 2000 #2000回投げて omega= 0.65 # 65%が表
genPriorA = omega * (kappa-2) + 1 #β分布のaパラメタに変換 genPriorB = ( 1.0 - omega ) * (kappa-2) + 1 # β分布のbパラメタに変換 genTheta<-rbeta(1,genPriorA,genPriorB) #パラメタ値生成
sampleN=74 # サンプルサイズを指定 sampleZ<-rbinom(1,size=sampleN,prob=genTheta) # データ生成 simulatedData<-c(rep(1,sampleZ),rep(0,sampleN-sampleZ)) # 0,1ベクトルに変換 simulatedData
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [71] 0 0 0 0
model_strings<-'
data{
int N;
int y[N];
}
parameters{
real<lower=0,upper=1> theta;
}
model{
theta ~ beta(1,1); #懐疑的な事前分布
y ~ bernoulli(theta);
}
'
library(rstan) model<-rstan::stan_model(model_code=model_strings) fit<-rstan::sampling(model,data=list(N=sampleN,y=simulatedData))
95%HDI=[0.492, 0.712]、ROPE=[0.48, 0.52]
ROPEの上限が95%HDIの下限を含むのでアウト!!!
Null_val<-0.5 # 空値 ROPE<- c(Null_val-0.02,Null_val+0.02) # ROPE theta.fit<-rstan::extract(fit)$theta # 事後分布のtheta値を抽出 thetaHDI_95<-quantile(theta.fit, probs=c(0.025,0.975)) #事後分布の95%HDIを算出 thetaHDI_95
## 2.5% 97.5% ## 0.4915092 0.7118824
95%HDIの幅=0.220, 望ましい幅(精度)=0.2なのでアウト
(thetaHDI_95[2]-thetaHDI_95[1])<0.2
## 97.5% ## FALSE
空値外しも精度もいずれの目標も事後分布から求められる
あとはこのフローを繰り返す(1000回とか任意の回数反復)
何回目標達成したかをカウントし目標達成確率を求める=検定力
# コイントスのstanコード
model_strings<-'
data{
int N;
int y[N];
}
parameters{
real<lower=0,upper=1> theta;
}
model{
theta ~ beta(1,1);
y ~ bernoulli(theta);
}
'
# stanコードをコンパイル
library(rstan)
model<-rstan::stan_model(model_code=model_strings)
# データ生成関数
simData.binom<-function(omega, kappa,sampleN){
#kappa= 2000 #2000回投げて
#omega= 0.65 # 65%が表
#sampleN=74 # サンプルサイズを指定
genPriorA = omega * (kappa-2) + 1 #β分布のaパラメタに変換
genPriorB = ( 1.0 - omega ) * (kappa-2) + 1 # β分布のbパラメタに変換
genTheta<-rbeta(1,genPriorA,genPriorB) #パラメタ値生成
sampleZ<-rbinom(1,size=sampleN,prob=genTheta) # データ生成
simulatedData<-c(rep(1,sampleZ),rep(0,sampleN-sampleZ)) # 0,1ベクトルに変換
return(simulatedData)
}
#目標達成関数
goalAchievedForSample <- function(fit,ROPE,W) {
thetaROPE = ROPE
thetaHDImaxWid= W
# 95%HDIを算出
thetaHDI<-as.vector(quantile(rstan::extract(fit)$theta,probs=c(0.025,0.975)))
# 目標: 空値ROPE除外
goalAchieved =list()
goalAchieved = c(goalAchieved,
"ExcludeROPE"=(thetaHDI[1]>thetaROPE[2]|thetaHDI[2]<thetaROPE[1]))
# 目標: 精度
goalAchieved = c(goalAchieved,
"NarrowHDI"=(thetaHDI[2]-thetaHDI[1]<thetaHDImaxWid))
return(goalAchieved)
}
fits<-list()
nSimulatedDataSets<-10 #シミュレーション数を指定
goalTally=NULL
for (simIdx in 1:nSimulatedDataSets){
simulatedData<-simData.binom(0.70,2000,74) #仮説パラメタ分布のモード (ω)と集中度 (κ), サンプルサイズ
standata<-list(N=length(simulatedData), y=simulatedData)
options(mc.cores = parallel::detectCores())
fit<-rstan::sampling(model,data=standata, cores = getOption("mc.cores", 4L))
fits[[simIdx]]<-fit #シミュレーション数を多くするときには外す
goalAchieved<-goalAchievedForSample(fit,c(0.48,0.52),0.2) #ROPE, widthを指定
if(!exists("goalTally")) {
goalTally=matrix(nrow=0, ncol=length(goalAchieved))
}
goalTally <- rbind(goalTally, goalAchieved)
}
for (goalIdx in 1:NCOL(goalTally)){
goalName=colnames(goalTally)[goalIdx]
goalHits=sum(unlist(goalTally[, goalIdx]))
goalAttempts=nrow(goalTally)
goalEst=goalHits/goalAttempts
show(paste0(goalName,": est.power=",round(goalEst,3)))
}
少ない回数で確認をしながら
固まってきたら回数増やす
joyplot使って見たかっただけ
## [1] "ExcludeROPE: est.power=1" ## [1] "NarrowHDI: est.power=0.4"
## [1] "ExcludeROPE: est.power=0.89" ## [1] "NarrowHDI: est.power=0.39"
反復回数が少ない間は、stanオブジェクトを全保存しても良いけど、反復が多くなると爆弾が出てくるのが怖いので、いちいち保存しない、目標達成したかどうかだけ保存していく
100回で少ない回数だけど、教科書と概ね一致 (教科書は1000回)
シミュレーションは論文をみると数千〜数万回当たり前
| theta.0.6 | theta.0.65 | theta.0.7 | theta.0.75 | theta.0.8 | theta.0.85 | |
|---|---|---|---|---|---|---|
| power.0.7 | 238 | 83 | 40 | 25 | 16 | 7 |
| power.0.8 | 309 | 109 | 52 | 30 | 19 | 14 |
| power.0.9 | 430 | 150 | 74 | 43 | 27 | 16 |
| theta.0.6 | theta.0.65 | theta.0.7 | theta.0.75 | theta.0.8 | theta.0.85 | |
|---|---|---|---|---|---|---|
| power.0.7 | 91 | 90 | 88 | 86 | 81 | 75 |
| power.0.8 | 92 | 92 | 91 | 90 | 87 | 82 |
| power.0.9 | 93 | 93 | 93 | 92 | 91 | 89 |
ただし、精度の例数設計は、幅を少し狭めると必要な例数がかなり増える 幅をどれくらいに見積もるかは、実質科学的観点から
誰でもエネルギーフィールド(EF)なるものがある (ん?)
病気はそれがその患部に凝固することによる (ん?)
直接触れない範囲で手を患者に近づけ凝固したEFをほぐす (ん?)
タッチセラピストは、エネルギーフィールドを感知でき (ん?)
患者の皮膚に手を近づけると、何かを感じとれる (調べます)
\(\omega\): 正答率の集団平均
\(\omega〜beta(A, B)\)
\(\kappa-2\): 集中度の集団平均
\(\kappa-2〜gamma(S, R)\)
\(\theta_s\): 個人の正答率
\(\theta_s〜beta(\omega(\kappa-2)+1, (1-\omega)(\kappa-2)+1)\)
\(\ y_{i|s}\):個人の個々の試行の回答
\(\ y_{i|s}〜berunoulli(\theta_s)\)
実データを使ったベイズ推定は、9章参照
理想的なデータサンプルでベイズ推定し、モデルパラメタの仮説的な事後分布を得る
その事後分布のパラメタを使ってデータセットを作って、目標達成確率を求めたい
各パラメタの不確実性を考慮した設計
# kruschkeの便利関数を使う為に読み込んどく
rm(list=ls())
source("script/Stan-Ydich-XnomSsubj-MbinomBetaOmegaKappa.R")
# 理想的な仮説を設定:
idealGroupMean = 0.65 # 群の平均回答率
idealGroupSD = 0.07 # 回答率の標準偏差
idealNsubj = 100 # 対象者 対象者を多くすれば、仮説の確信度は高くなる
idealNtrlPerSubj = 100 # 試行数 試行数を多くすれば、仮説の確信度は高くなる
# 理想的な仮説からパラメタ値をランダムに生成 betaAB = betaABfromMeanSD( idealGroupMean , idealGroupSD ) #平均とSDをベータ分布のパラメータに変換 theta = rbeta( idealNsubj , betaAB$a , betaAB$b ) # 仮説のパラメータ分布をrbetaで発生 # 理想とする平均、標準偏差にマッチするように少し補正: theta = ((theta-mean(theta))/sd(theta))*idealGroupSD + idealGroupMean theta[ theta >= 0.999 ] = 0.999 # must be between 0 and 1 theta[ theta <= 0.001 ] = 0.001 # must be between 0 and 1 # 仮説のパラメタ値に近似するデータを生成 z = round( theta*idealNtrlPerSubj )
# 1人100回の回答を100名分, y=回答, sはID. theta値に基づいて、0,1データに変換
dataMat=matrix(0,ncol=2,nrow=0,dimnames=list(NULL,c("y","s")))
for ( sIdx in 1:idealNsubj ) {
yVec = c(rep(1,z[sIdx]),rep(0,idealNtrlPerSubj-z[sIdx]))
dataMat = rbind( dataMat , cbind( yVec , rep(sIdx,idealNtrlPerSubj) ) )
}
idealDatFrm = data.frame(dataMat)
head(idealDatFrm$y,200)
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 ## [71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 ## [106] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## [141] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## [176] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
head(idealDatFrm$s,200)
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ## [71] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 ## [106] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ## [141] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ## [176] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
y = as.numeric(idealDatFrm[,"y"])
s = as.numeric(idealDatFrm[,"s"]) # ensures consecutive integer levels
z = aggregate( y , by=list(s) , FUN=sum )$x #個人の正答数カウント
N = aggregate( rep(1,length(y)) , by=list(s) , FUN=sum )$x #個人の試行数
Nsubj = length(unique(s)) # 被験者数
# Stanに渡すデータセット
dataList = list(
z = z ,
N = N ,
Nsubj = Nsubj
)
modelString = "
data {
int<lower=1> Nsubj ;
int<lower=0> z[Nsubj] ;
int<lower=0> N[Nsubj] ;
}
parameters {
real<lower=0,upper=1> theta[Nsubj] ; // individual prob correct
real<lower=0,upper=1> omega ; // group mode
real<lower=0> kappaMinusTwo ; // group concentration minus two
}
transformed parameters {
real<lower=0> kappa ;
kappa <- kappaMinusTwo + 2 ;
}
model {
omega ~ beta( 1 , 1 ) ;
kappaMinusTwo ~ gamma( 0.01 , 0.01 ) ; // mean=1 , sd=10 (generic vague)
// kappaMinusTwo ~ gamma( 1.105125 , 0.1051249 ) ; # mode=1 , sd=10
theta ~ beta( omega*(kappa-2)+1 , (1-omega)*(kappa-2)+1 ) ; // vectorized
for ( s in 1:Nsubj ) {
z[s] ~ binomial( N[s], theta[s] ) ;
}
}
"
model<-rstan::stan_model(model_code=modelString)
#initListはKruschkeの自作関数 (初期値をいい感じにふってくれる)
parameters = c( "theta","omega","kappa") # The parameters to be monitored
burnInSteps = 500 # Number of steps to burn-in the chains
nChains = 4 # nChains should be 2 or more for diagnostics
numSavedSteps =2000
thinSteps = 20
# Translate to C++ and compile to DSO:
#stanDso <- stan_model( model_code=modelString )
# Get MC sample of posterior:
stanFit <- rstan::sampling( object=model ,
data = dataList ,
pars = parameters , # optional
chains = nChains ,
iter = ( ceiling(numSavedSteps/nChains)*thinSteps
+burnInSteps ) ,
warmup = burnInSteps ,
thin = thinSteps)
#ステップ4: パラメタ値からランダムにデータセットを生成
#ステップ5: 生成されたデータセットでMCMC->目標達成の有無をカウント
# 一つのデータセットのMCMCで目標が達成されたかを確認する関数
goalAchievedForSample = function( data ) {
nullROPE = c(0.48,0.52)
HDImaxWid = 0.2
# MCMCをstanで実行, Kruschkeの自作関数
mcmcCoda = genMCMC( data=data , saveName=NULL ,
numSavedSteps=5000 , thinSteps=2 )
mcmcMat = as.matrix(mcmcCoda)
# HDIの算出:
HDImat = apply( mcmcMat , 2 , "HDIofMCMC" )
show( HDImat[,1:5] )
# 目標達成をレコードする箱:
goalAchieved = list()
# 目標: omegaのROPE外し:
goalAchieved = c( goalAchieved ,
"omegaAboveROPE"=unname( HDImat[1,"omega"] > nullROPE[2] ) )
# 目標: omega の精度固め
goalAchieved = c( goalAchieved ,
"omegaNarrowHDI"=unname( HDImat[2,"omega"]-HDImat[1,"omega"]
< HDImaxWid ) )
# 目標: 少なくとも1つのthetaがROPEのより大きく、ROPEより小さいものはなし
thetaCols = grep("theta",colnames(HDImat)) # column indices of thetas
goalAchieved = c( goalAchieved ,
"thetasAboveROPE"= (any(HDImat[1,thetaCols] > nullROPE[2])
& !any(HDImat[2,thetaCols] < nullROPE[1])))
# 目標: 全てのthetaの精度固め
goalAchieved = c( goalAchieved ,
"thetasNarrowHDI"= all( HDImat[2,thetaCols]
- HDImat[1,thetaCols]
< HDImaxWid ) )
# More goals can be inserted here if wanted...
# Return list of goal results:
return(goalAchieved)
}
stanDso <- model
# シミュレーションの設定
#Nsubj = 2*7 ; NtrlPerSubj = 47 # 個人の数を大きく
Nsubj = 7 ; NtrlPerSubj = 2*47 # 試行数を大きく
# シミュレーション回数を設定:
nSimulatedDataSets = min(100,NROW(mcmcMat))
# Run the simulated experiments:
simCount=0
if (exists("goalTally")) rm(goalTally)
for ( simIdx in ceiling(seq(1,NROW(mcmcMat),length=nSimulatedDataSets)) ) {
simCount=simCount+1
cat( "\n\n==================== Simulation",simCount,"of",nSimulatedDataSets,
"====================\n\n" )
# 理想的なデータから生成したomega,kappaから値を一つ:
genOmega = mcmcMat[simIdx,"omega"]
genKappa = mcmcMat[simIdx,"kappa"]
# そこからtheta値を生成:
genTheta = rbeta( Nsubj , genOmega*(genKappa-2)+1 , (1-genOmega)*(genKappa-2)+1 )
# そこからデータセットを生成:
dataMat=matrix(0,ncol=2,nrow=0,dimnames=list(NULL,c("y","s")))
for ( sIdx in 1:Nsubj ) {
z = rbinom( 1 , size=NtrlPerSubj , prob=genTheta[sIdx] )
yVec = c(rep(1,z),rep(0,NtrlPerSubj-z))
dataMat = rbind( dataMat , cbind( yVec , rep(sIdx,NtrlPerSubj) ) )
}
# MCMC実行、目標達成カウント:
goalAchieved = goalAchievedForSample( data.frame(dataMat) )
# Tally the results:
if (!exists("goalTally")) { # if goalTally does not exist, create it
goalTally=matrix( nrow=0 , ncol=length(goalAchieved) )
}
goalTally = rbind( goalTally , goalAchieved )
}
load("data/Tally2.Rdata")
# 目標達成のカウント
for ( goalIdx in 1:NCOL(goalTally) ) {
# Extract the goal name for subsequent display:
goalName = colnames(goalTally)[goalIdx]
# Compute number of successes:
goalHits = sum(unlist(goalTally[,goalIdx]))
# Compute number of attempts:
goalAttempts = NROW(goalTally)
# Compute proportion of successes:
goalEst = goalHits/goalAttempts
# Compute HDI around proportion:
goalEstHDI = HDIofICDF( qbeta ,
shape1=1+goalHits ,
shape2=1+goalAttempts-goalHits )
# Display the result:
show( paste0( goalName,
": Est.Power=" , round(goalEst,3) ,
"; Low Bound=" , round(goalEstHDI[1],3) ,
"; High Bound=" , round(goalEstHDI[2],3) ) )
}
# シミュレーションの設定 #Nsubj = 2*7 ; NtrlPerSubj = 47 # 個人の数を大きく Nsubj = 7 ; NtrlPerSubj = 2*47 # 試行数を大きく
## [1] "omegaAboveROPE: Est.Power=0.37; Low Bound=0.28; High Bound=0.466" ## [1] "omegaNarrowHDI: Est.Power=0.31; Low Bound=0.226; High Bound=0.404" ## [1] "thetasAboveROPE: Est.Power=1; Low Bound=0.971; High Bound=1" ## [1] "thetasNarrowHDI: Est.Power=0.87; Low Bound=0.795; High Bound=0.926"
集団レベルの推定の精度悪く, 個人レベルの推定精度良し
# シミュレーションの設定 Nsubj = 2*7 ; NtrlPerSubj = 47 # 個人の数を大きく #Nsubj = 7 ; NtrlPerSubj = 2*47 # 試行数を大きく
## [1] "omegaAboveROPE: Est.Power=0.98; Low Bound=0.938; High Bound=0.997" ## [1] "omegaNarrowHDI: Est.Power=0.97; Low Bound=0.923; High Bound=0.993" ## [1] "thetasAboveROPE: Est.Power=1; Low Bound=0.971; High Bound=1" ## [1] "thetasNarrowHDI: Est.Power=0.28; Low Bound=0.199; High Bound=0.372"
集団レベルの推定の精度良く, 個人レベルの推定精度悪い
> 有意差出なかったら論文出せないよ〜
> 空値についてなんと言われようと真実の愛に突き動かされている
ガンジー先輩愛を語る
' 愛に基づく力は、恐怖に基づく力よりも1000倍有効で永続するんだぜ '
θ=0.50: 初期のシーケンスで空値を誤って棄却
θ=0.65: 初期のシーケンスで誤って帰無仮説を採択
θ=0.50: 初期のシーケンスで空値を誤って棄却
θ=0.65: 初期のシーケンスで誤って帰無仮説を採択
θ=0.50: 初期のシーケンスで帰無仮説を棄却しない
θ=0.65: 初期のシーケンスで誤って帰無仮説を採択しない 初期の不確実性が適切に反映されている
ROPEの80%になったら中止 ROPEを0.45~0.55に設定するとHDIの臨界幅は0.08 0.08以下になるまでデータ収集が続く
そのポイントでHDIとROPEを比較し、帰無仮説を採択するか否か判断可能
帰無仮説が真であっても、常に帰無仮説を棄却することになりθを過剰推定しがち
帰無仮説が真である場合はフォールスアラームを100%防げるが、帰無仮説が偽である場合には誤って帰無仮説を採択しがち θが極端になりがち
フォールスアラームを100%防げて、帰無仮説が偽である場合に誤って帰無仮説を採択することもない
サンプルサイズ大きくなりがち θが極端になることある
フォールスアラームを100%防げて、帰無仮説が偽である場合に誤って帰無仮説を採択することもない
判断しないことありがち サンプルサイズ大きくなりがち θが極端になることはない(パラメタの値に測定の精度が依存しない限り)
library(BEST)
library(rjags)
# 1. Generate idealised data set:
proData <- BEST::makeData(mu1=20, sd1=5, mu2=10, sd2=5, nPerGrp=1000,
pcntOut=10, sdOutMult=2.0, rnd.seed=123)
# 2. Generate credible parameter values from the idealised data:
proMCMC <- BESTmcmc(proData$y1, proData$y2, numSavedSteps=2000)
N1plan <- N2plan <- 20
powerPro <- BESTpower(proMCMC, N1=N1plan, N2=N2plan,
ROPEm=c(-1.5,1.5), ROPEsd=c(-2,2), ROPEeff=c(-0.5,0.5),
maxHDIWm=15.0, maxHDIWsd=10.0, maxHDIWeff=1.0, nRep=10,
showFirstNrep=5,saveName="test")